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ABSTRACT 

We present mid-infrared spectral maps of the NGC 1333 star forming region, 
obtained with the the Infrared Spectrometer on board the Spitzer Space Telescope. 
Eight pure H2 rotational lines, from 5(0) to 5(7), are detected and mapped. The 
H2 emission appears to be associated with the warm gas shocked by the multiple 
outflows present in the region. A comparison between the observed intensities 
and the predictions of detailed shock models indicates that the emission arises 
in both slow (12 - 24 km/s) and fast (36 - 53 km/s) C-type shocks with an 
initial ortho-to-para ratio < 1. The present H2 ortho-to-para ratio exhibits a 
large degree of spatial variations. In the post-shocked gas, it is usually about 2, 
i.e. close to the equilibrium value (~ 3). However, around at least two outflows, 
we observe a region with a much lower (~ 0.5) ortho-to-para ratio. This region 
probably corresponds to gas which has been heated-up recently by the passage 
of a shock front, but whose ortho-to-para has not reached equilibrium yet. This, 
together with the low initial ortho-to-para ratio needed to reproduce the observed 
emission, provide strong evidence that H2 is mostly in para form in cold molecular 
clouds. The H2 hues are found to contribute to 25 - 50% of the total outflow 
luminosity, and thus can be used to ascertain the importance of star formation 
feedback on the natal cloud. From these lines, we determine the outflow mass 
loss rate and, indirectly, the stellar infall rate, the outflow momentum and the 
kinetic energy injected into the cloud over the embedded phase. The latter is 
found to exceed the binding energy of individual cores, suggesting that outflows 
could be the main mechanism for core disruption. 

Subject headings: astrochemistry — stars: formation — ISM: abundances — 
ISM: molecules — ISM: individual () 
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Introduction 



The early phases of the formation of a star are characterized by the presence of bipolar 
outflows, which, as the material falls onto the young star, eject matter at high speeds (several 



tens of km/s) to distances up to tens of parsec away (e.g. Arce et al. 2007). These outflows 



have a great impact on the surrounding gas. They can clear out cavities in their natal 
molecular clouds, and may eventually cause cloud disruption. As they propagate through 
molecular clouds, they create shock waves that compress, heat up and alter the chemical 



composition of the gas. Theoretical studies of shock waves in the ISM (see Draine &; McKee 



1993 , for a review) predict that their nature depends on the shock velocity and the intensity 
of the magnetic field. Slow shocks (< 45 km/s for typical values of the magnetic field) are 
of type C: the density, temperature and other quantities vary continuously through their 
passage. Furthermore, they are non-dissociative: molecules generally survive their passage, 
although the chemical composition of the gas can be significantly altered. On the other 
hand, faster shocks produce discontinuities in both the density and temperature of the gas, 
and may dissociate molecules on their passage. 

Observations of H2 pure rotational lines in shock regions can provide important con- 
straints on the physical and chemical conditions that prevail in these regions. For example, 
the H2 ortho-to-para ratio (hereinafter opr) is predicted to be greatly affected by the passage 
of a C-type shock. In the cold (pre-shock) gas, H2 is expected to be mostly in para form 



(Flower et al. 2006 Maret & Bergin 2007). In the shocked gas, reactive collisions with H 



atoms can convert P-H2 into 0-H2 ( Timmermann|[l998 Wilgenbus et al.||2000 ). This reaction 
has a relatively high activation barrier (~ 4000 K), and therefore the efficiency of the conver- 



sion depends critically on the temperature that is reached in the shock. Indeed, Wilgenbus 



et al. (2000 ) showed that effective para to ortho H2 conversion takes place in the temperature 



interval 700 — 1300 K, while Kristensen et al. (2007) refined this result in to the interval 
800 — 3200 K. In addition, the efficiency of the conversion depends on the abundance of H 
atoms in the gas. Fast shocks produce higher H abundances, which makes the conversion 
faster. 

Molecular hydrogen rotational lines have been observed in several star forming regions 
with the Infrared Space Observatory and, more recently, with the Spitzer Space Telescope 
(Werner et al. 2004). Using the Short Wavelength Spectrometer on board ISO, Neufeld 
et al. (1998) observed five pure rotational lines (from — 5(1) to — 0S'(5), hereafter 
simply referred as S{1) and S{5), respectively) towards the Herbig-Haro 54 outfiow (HH 54). 
Using these multiple ortho and para lines, they determined simultaneously the rotational 
temperature and the ortho-to-para ratio. Interestingly, they measured an ortho-to-para 
ratio of 1.3, which is significantly lower than the equilibrium value (~ 3) at the observed 
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gas temperature (~ 650 K). Lefloch et al. (2003) used the ISOCAM camera to map the 



H2 S(2) to S{y) towards the HH 2 outflow, and found important spatial variations in the 



opr (between 1.2 and 2.5) around the HH 2 object. Neufeld et al. (2006) used the Infrared 



Spectrometer (IRS, Houck et al. 2004) on board Spitzer to map the opr towards HH 54 and 
HH 7-11, and also observed a low, non-equilibrium opr, with important spatial variations. 
The fact that the opr is out-of-equilibrium suggests that the ortho-to-para conversion is 
relatively slow in these regions: the observed opr was interpreted as a fossil of an earlier 
epoch when the gas was cooler. 

This paper presents complete 5.2 — 36.5 fim spectral maps obtained with the IRS towards 
the NGC 1333 region in the Perseus cloud. The observations cover a region of roughly 6' x 10', 
i.e. 0.4 pc X 0.5 pc assuming a distance of 220 pc ( Cernis||1990 ). To our knowledge, these are 
the largest complete spectral maps obtained with the IRS. NGC 1333 contains many YSOs, 
revealed by sub-millimeter ( Sandell fc Knee||2001 ), millimeter (Lefloch et al. 1998) and mid- 



infrared continuum maps (Gutermuth et al. 2008). It also contains numerous molecular 



outflows, detected for example in CO line emission (Knee & Sandell 2000). Here we focus 



on H2 pure rotational line emission, from S{0) to S{7). These lines are used to map the 
H2 rotational temperature and opr over a large region, thus permitting detailed study of 
the impact of outflows on the surrounding cloud. The paper is organized as follow. In ^ 
we present the observations and discuss the data reduction. The spectra and maps that 
we obtained are presented in ^ In ^ we describe the analysis procedure that we used to 
derive the H2 column density, rotational temperature and opr from our observations. In ^ 
we compare our observations with the predictions of shock models, we compare the variation 
of the opr with the temperature, and we discuss the impact of the outflows on the cloud. 
Our conclusions are presented in ^ 



2. Observations 

NGC 1333 was observed using the IRS during Cycle 2 of the General Observer program, 
in March and September 2006. The Long-High (LH), Short-High (SH) and Short-Low (SL) 
modules of the IRS were used, providing a complete spectral coverage from 5.2 to 36.5 /xm. 
The spectral resolution of these observations {R = A/ A A) ranges from 64 to 128 for the SL 
mode, and is ~600 for the high resolution modules (SH and LH). The half-power beam size 



of the instrument ranges from 3" at 5.2 /im to 10" at 38 fxia (Neufeld et al. 2006) 



The SH observations consist of nine different Astronomical Observations Requests (AORs), 
each covering a rectangular region of ~ 142" x 120". Each of these AORs was obtained by 
moving the IRS slit in both perpendicular and parallel directions of the slit length, with a 
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half-slit width stepping in perpendicular direction, and full-slit length in the parallel direc- 
tion. Our SH observations do not cover the tip of the HH 7-11 outflow; therefore we have 
used archival data of this region obtained in February 2004 as part of the IRAC Guaran- 
teed Time program, and that consists of two AORs of 58" x 76" each. Details on these 



data can be found in Neufeld et al. (2006). When merged, the SH AORs cover a region of 
~ 6' X 10', roughly centered on SVS 13 (a = 17^22'"38.2'^ and S = -23°49'34.0"; J2000). 
The LH observations consist of six AORs of ~ 200" x 273" each. As for the SH observations, 
these were obtained with a half-slit width stepping in perpendicular direction, and full-slit 
length in the parallel direction. These observations cover a region of ~ 8' x 9'. Finally, the 
SL observations consist of six AORs of ~ 202" x 222", each obtained with a full-slit width 
stepping in the perpendicular direction, and half-slit width stepping in the parallel direction. 
These observations cover a region of ~ 6' x 14'. 

All the data, including the archival ones, have been reduced using the latest version 
(17.2) of the Spitzer Science Center pipeline. Further reduction was performed using the 
SMART software package ( Higdon et al.||200^ ), supplemented by the IDL routines described 
in Neufeld et al. (2006). These routines remove the bad pixels in the LH and SH observations. 



extract spectra for each pixel along the slit of the different modules, spatially re-sample 
the spectra on a regular grid (with a spatial resolution corresponding to the pixel spacing, 
i.e. 2.3", 1.8" and 4.5" for the SH, SL and LH modes, respectively), and finally produce 
continuum-subtracted spectral maps by fitting a Gaussian for each spectral line. A first 
look at our LH and SH spectral maps revealed the presence of stripes in the direction 



perpendicular to the slit length. As discussed in Neufeld et al. (2007), these stripes are due 



to the bad pixel extraction routine, which interpolates missing pixel values in the dispersion 
direction. Therefore, if a bad pixel is close to the central wavelength of a given line, the 
intensity measured at a given position along the slit is consistently underestimated, creating 



the observed stripes. Following Neufeld et al. (2007), we have corrected this effect by applying 



a correction factor for each line and each position along the slit (see Neufeld et al. 2007 for 



details on how this correction factor is determined). This technique was found to remove 
the stripes from our maps quite efficiently. Note that our LH observations were not affected 
by this problem because they were obtained with a half-slit length stepping in the parallel 
direction, as opposed to the SH and SL data which were obtained with a full-slit length 
stepping in that direction. Therefore, in LH mode, any given position on the sky was 
observed twice, and potential stripes were washed-out when averaging the spectra together. 
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3. Results 
3.1. Spectra 

Fig. [T] to [3] present average spectra obtained on several positions of our maps. These 
spectra were obtained by averaging all the spectra observed within a 15" FWHM Gaussian 
aperture around a given position, in order to increase the signal-to-noise ratio. We detect 
eight pure rotational lines - from S{0) to S{7) - towards HH 7. On the other hand, no H2 
emission is detected towards the NGC 1333-IRAS 4A and 4B Class protostars. Several 
fine structure atomic lines are also detected, such as the Fe II ^-^7/2 — ^ -F9/2 (17.9 /im), the 
S I -3 P2 (25.2 fxm), Fe II ^^7/2 P'9/2 (26.0 /zm) and Si II ^po^^ P?/2 (34-8 fxm) lines. 
Table [T] gives the intensities measured towards several positions for these lines. Thanks to 
an overlap in the wavelength range covered by the SL and SH modules, the H2 S{2) intensity 
was measured with each module, and was found to agree within ~ 20%. 

Several PAH features - at 6.2, 7.2, 8.6, 11.2 and 12.8 fim - as well as the 9.7 /xm silicates 
absorption band are clearly detected towards IRAS 8. Finally, the CO2 ice feature at 15.2 fim 
is seen on IRAS 2. A discussion of all the detected spectral features is beyond the scope of 
this paper; in the following, we focus on the interpretation of the H2 line emission, and we 
postpone the discussion of other lines and spectral features to forthcoming papers. 



3.2. Maps 

Fig. 4 a-i show maps of the H2 5(0) to 5(7) lines. H2 5(1), 5(2), 5(3), 5(4) and 5(5) 
emission is readily detected along several outflows. The H2 5(0), 5(6) and 5(7) maps have 
a lower signal-to-noise ratio, and emission of these lines is barely detected at this spatial 
resolution (4.5 and 1.8" for the LH and SL modes, respectively). In the following, we discuss 
the morphology of the H2 emission along the different outflows. 



3.2.1. SVS 13 and HH 7-11 



The Class I source SVS 13 ( Strom et al.||1976" ) is associated with a chain of Herbig-Haro 
objects, HH 7-11. These objects are excited by a high-velocity molecular outflow, which is 
detected in CO and SiO rotational lines ( [Bachiller et"aL]|1998a[ [2000} [Knee fc Sandell|[2000| ) , 
as well as the H2 1 — 5(1) ro- vibrational line at 2.12 /im (Aspin et al. 1994 Hodapp & 



Ladd 1995 Khanzadyan et al. 2003). A second outflow, roughly orientated along a north- 
south axis and driven by the SVS 13B Class source, is detected in SiO J = 2 — 1 emission 
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(iBachiller et al.||1998bl) 



The south-east lobe of the HH 7-11 outflow is clearly detected in our S{1) map (Fig. 
4 a). This map does not cover the north-west lobe, but emission from this lobe is seen in the 



S'(2) and S{3) maps (Fig. 4 b and 4 c). H2 S{1) and S{3) emission is also detected along the 
north-south outflow originating from SVS 13B. The position angle (from the north to east) 
of the H2 emission is about 160°. The emission extends as far as ~ 3' south of SVS 13B, up 
to the tip of the NGC 1333-IRAS 4A outflow (see below). Some weak emission along this 
flow is also seen in the H2 >S'(2) line. Diffuse S{1) emission extends to the north along the 



same axis, and may be associated with this outflow as well (see ^3.2.5) 



3.2.2. NGC 1333-IRAS 2 



NGC 1333-IRAS 2 ( [Jennings et al.|1987| hereafter IRAS 2) is a binary system composed 
of the IRAS 2A and IRAS 2B Class protostars. This system is associated with two outflows 
that both originate ~ 6" west of IRAS 2 A ( Knee & Sandell 2000 ) . The first one is roughly 



orientated along a north-south axis, and is detected in CO (Liseau et al. 1988 Sandell et al. 



1994| and H2 1 - 5(1) ( [Hodapp fc Ladd|[l995| . The second outflow is roughly oriented 
along an east-west axis (position angle ~ 104°) and is detected in CO, CH3OH and SiO 



(Sandell et al. 1994 Bachiller et al. 1998a; Knee & Sandell 2000; J0rgensen et al. 2004). 



In our maps, we observe a peak of 5'(1), 5'(2) and S{3) emission east of IRAS 2B along 
the axis of the east-west outflow. The position of this peak corresponds to the peak of the CS 
and CH3OH emission, and presumably arises from a bow-shock ( Langer et al.[[1996 Bachiller 



et al. 1998a). No emission is seen west of the source, but our maps have limited coverage in 



that direction. Two II2 emission spots are also seen north and south of IRAS 2A. These two 
peaks lie at a position angle of ~ 12°, consistent with the orientation of the H2 1 — >S'(1) 
jet, but slightly different from the position angle of the CO jet (~ 25° Sandell et al. 1994). 
No emission is seen towards the protostars. 



3.2.3. NGC 1333-IRAS 4 



NGC 1333-IRAS 4 is also a binary system, composed of two Class protostars, IRAS 4A 
and IRAS 4B. IRAS 4A drives a powerful and well coUimated outflow that extends roughly 
over 4', and is detected in H2 1 — S{1) ( [Hodapp &: Ladd 1995) as well as CO and SiO lines 



(Blake 1996 Knee &; Sandell 2000). A more compact outflow, with an inclination close to 



90° (i.e. almost perpendicular to the plane of the sky), is driven by IRAS 4B (Blake 1996). 
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We detect H2 S{1) emission at the southwestern tip of the IRAS 4A outflow. Unfortu- 
nately, our SH maps do not cover other parts of the flow. However, S{2) and S{3) emission 
peaks are seen at to the north-east and south-west of IRAS 4A. No H2 is detected from the 
IRAS 4B outflow, nor from IRAS 4A and IRAS 4B protostars themselves. 



3.2.4. NGC 1333-IRAS 7 and HH 6 



NGC 1333-IRAS 7 is associated with two sub-millimeter sources, SM 1 and SM 2 



(Sandell & Knee 2001). The former is coincident with the cm source VLA 27 (Rodriguez 



et al.|1999 ), as well as with a water maser ( Henkel et al.|1986 ). A molecular outflow, oriented 



along a east-west axis and probably driven by SM 1, is detected in CO (Liseau et al. 1988 



Knee fc Sandell||2000 ). This outflow is the driving source for the chain of HH objects located 
east of SM 1, which include HH 6. Our maps reveal bipolar H2 emission orientated along 
the same axis as the CO outflow. 



3.2.5. NGC 1333-IRAS 6 and HH 12 



NGC 1333-IRAS 6 is coincident with a PMS star (SVS 107) and the continuum cm 



source VLA 42 (Rodriguez et al. 1999). An elongated structure is also seen towards that 



source in the sub-mm maps of Sandell & Knee (2001). west of NGC 1333-IRAS 6 is a chain 



of HH objects that include HH 12. Knee & Sandell (2000) suggested that these objects are 



excited by the outflow originating from SVS 13B. 

In our H2 maps we see some diffuse emission west of IRAS 6 that extends from the 
south down to the north west lobe of the SVS 13 outflow. Although this diffuse emission 
appears broadly aligned with the axis of the SVS 13B, it is difficult to say from the present 
data if the H2 emission arises in the shocked gas from this outflow. 



3.2.6. NGC 1333-IRAS 8 



NGC 1333-IRAS 8 is part of the optical nebulosity north of the cloud (see Fig. 4 b). 
It is associated with a PMS star, SVS 3. Although the IRS spectrum is dominated by PAH 
features (see Fig. [T]), we also detect extended H2 5'(2), S'(3) and S{h) (but no 5'(4)) to the 
south east of the source. 
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4. Analysis 

4.1. Rotational diagrams 

In Fig. [sj we show rotational diagrams obtained for 15" FWHM Gaussian apertures 
centered at several positions along the outflows: HH 7, IRAS 4A-SW, IRAS 2A-S, IRAS 7 
and SVS 13B-S. These positions are indicated on the H2 maps (Fig. 4 a-i). Line intensities 
were corrected for extinction using the standard = 3.1) galactic dust opacity curve from 



Weingartner & Draine (2001), assuming = 2 for all sources (the value measured towards 
HH 8; |Gredel|1996j ). In each dia gram and for each line, we have plotted ln{Nu/ {guQs))-, where 
Nu is the column density in the upper level, Qu is the rotational degeneracy (equal to 2 J + 1) 
and Qs is the spin degeneracy (3 for ortho, and 1 for para transitions). If the ortho-to-para 
ratio is equal to its high temperature limit {i.e. 3), and if H2 transitions are thermalized at 
a single temperature, then the observed values would line up on a straight line. Instead, the 
diagram show a "zig-zag" pattern, indicating that the ortho-to-para ratio is lower than 3. In 
addition, the rotational diagrams show a positive curvature; the observed values cannot be 
flt with a single excitation temperature. This suggests that the observed H2 emission arises 
in a mixture of gas with different temperatures, as already observed in HH 7 and HH 54 



by Neufeld et al. (2006). Following these authors, we have fltted the observations with two 
gas components (herein after the "warm" and the "hot" components). Note that this two 
temperature model is an approximation; the gas has more likely a continuous temperature 
distribution. The total H2 column density, ortho-to-para ratio and excitation temperature in 
each gas component were kept as free parameters, and were adjusted to fit the observations 
using a minimization routine. These three parameter can be constrained quite independently 
from each other. The ortho-to-para ratio is constrained from the "zig-zag" pattern in the 
rotational diagram. The rotational temperature is determined from the slope of the rotational 
diagram. Finally, the total H2 column density is determined from the ordinate of the points 
in the diagram. This procedure produced good fits to the data; in particular, no clear 
departures from the LTE is seen, even for the highest energy lines. The critical densities - 
above which collisional excitation dominates over radiative de-excitation, and energy levels 
are populated according to Boltzmann's statistic - range between 5 and 4 x 10^ cm~^ for 



the 5'(0) and S{7) lines respectively (assuming a kinetic temperature of 1000 K; Le Bourlot 



et al.|1999 ). Therefore the highest energy lines (>S'(6) and S'(7)) must arise in relatively dense 



5 oT>-,-3^ 



gas (> 10 cm 

The best fit parameters for HH 7, IRAS 4A-SW, IRAS 2A-S, IRAS 7, SVS 13B-S 
are given in Table [2] The rotational temperature towards these sources is between ~ 300 
and ~ 600 K for the warm component, and ~ 1000 - 1500 K for the hot component. We 
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measure an opr of 0.3 - 0.7 and 1.4 - 2.1 in the warm and hot components, respectively]^ 
For both components, the opr is much lower than its equilibrium value for the measured 
(rotational) gas temperatures; for kinetic temperatures higher than 300 K, the equilibrium 
value is essentially 3. For all sources, we measure a higher opr in the hot component than in 
the warm component. Similar behavior was noted previously by Neufeld et al. (2006, 2007). 
The dependence of opr on temperature is discussed in further detail below. The H2 column 
densities are typically ~ 2 — 6 x 10^^ cm^ and ~ 3 — 6 x 10^^ cm^ in the warm and hot 
components, respectively. Thus, the column density of the hot component is more than an 
order of magnitude smaller than the column density of the warm component. 



4.2. H2 column density, rotational temperature and ortho-to-para ratio maps 

Thanks to their large spatial coverage, our observations essentially map the ortho-to- 
para ratio, rotational temperature, and column density of H2 in the core of the NGC 1333 
cloud. For this purpose, the H2 line maps were re-sampled on a common grid with a 5" 
spacing. Only regions of the sky that were observed with all three modules were considered. 
We then constructed a rotational diagram in each pixel of the resulting map. As mentioned 
above, the H2 rotational emission arises from at least two gas components. With three free 
parameters for each component (column density, opr and rotational temperature), we need at 
least six lines detected in a given pixel of the maps. In practice, few pixels in our maps meet 
this criterion. Therefore, when constructing the maps, only the H2 S'(O) to >S'(3) transitions 
were considered, which allows for greater spatial coverage of the region. Since the 5(0) to 
S'(3) line emission is dominated by the warm component (see Fig. ]5]), the column density, 
opr and rotational temperature derived using these lines correspond to that component. 

On Fig. ]6| ]7] and ]8] we provide maps of the H2 column density, rotational temperature 
and opr. These quantities could be measured along both lobes of the SVS 13 outflow and 
around IRAS 7. We could also map those at the south-west tip of the IRAS 4A outflow, as 
well at the northern, southern and western tips of the IRAS 2 outflows. Elsewhere in the 
map, we did not detect the H2 5'(0), 5(1), S{2) and 5(3) lines simultaneously, or the lines 
were detected with an insufficient signal-to-noise ratio. Where measured, some important 
variations in these quantities are seen. The H2 column density ranges from 3 x 10^^ to 10^° 
cm~^. The highest column densities are seen along the HH 7-11 outflow (see Fig. |6|). The 
rotational temperature ranges between 200 and 650 K, with a maximum reached close to the 



^The best fit parameters we obtain for HH 7 differ slightly from those reported by 'Neufeld et al. ('2006' 
on the same source. The reason is that we used new data in addition to those ofiNeufeld et al., (.2006,). 
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position of HH 7, as well as the IRAS 2N position (see Fig. [?]). Finally, the H2 ortho-to-para 
ratio ranges from 0.1 and 3 (Fig. Isl). 



Some important variations in both the rotational temperature and the opr are seen. 
Along the HH 7-11 outflow, between SVS 13 and HH 8, we measure a rotational temperature 
of ~ 400 K and an opr ranging between 1.5 and 2. Around HH 7, we measure a higher 
rotational temperature (~ 650 K), but a lower opr (~ 0.5). The presence of hot gas with a 



low opr in front of the HH 7-11 has already been observed by Neufeld et al. (2006). They 
suggested that this gas has been compressed and heated-up by the passage of the shock, but 
the ortho-to-para ratio has not reached LTE yet. Along the northwestern lobe of this flow, 
the spectra are noisier, and large variations from one pixel to the next. Around IRAS 4- 
SW, the rotational temperature and the opr were measured only in a few pixels making 
the comparison between these two quantities difficult. On the other hand, we observe some 
variation in the rotational temperature along the north and south lobes of the IRAS 2 
outflow. Like in HH 7-11, the rotational temperature is higher at the tip of this outflow. 
However, no corresponding variation is seen in the opr. Some relatively hot gas (~ 500 K) 
with a low opr (~ 1) is also observed around the eastern lobe of the IRAS 7 outflow. The 
variation of the opr as a function of the Trot are further discussed in Appendix lAl 



5. Model and Discussion 

5.1. Comparison of the observed emission with the predictions of a shock 

model 

In this section, we compare the observed H2 line emission towards several outflow po- 
sitions with the predictions of a stationary, planar C-type shock model that includes ortho 



para interconversion (Flower & Pineau des Forets 2003). The shock model has several free 



parameters: the pre-shock density n, the pre-shock velocity v, the initial (pre-shock) ortho- 
to-para ration opvi, and the intensity of the transverse magnetic fleld. The latter is expressed 

— 1/2 

as b = Bf Uq , where is the transverse magnetic fleld in jjG and uq is the pre-shock den- 
sity in cm~^. The model has been run for a grid of parameters (Kristensen et al. in prep), 
and the results of these computations are available on the wetj^ 

In Fig. |9| we compare the column densities obtained from our observations on aperture 
averaged positions with those predicted by the shock model. As a flrst approach, we have 
tried to reproduce the observed column densities with a single component shock. The value 



http : //www . strw . leidenuniv . nl/$\sim$kristensen/models 
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of h was set to 1, and the other parameters of the model were adjusted until a good agreement 
with the observations was obtained. We found that the observed 5(0), 5(1), 5(2) emission 
were fairly well fitted for a pre-shock density of 10^ cm~^, an initial o'pr of 0.01 or 1, a 
pre-shock velocity of 12-24 km/s, and a beam filling factor // of typically 10%. However, 
this model greatly underestimates the observed lines with higher energies. In other words, 
the temperature reached in a slow shock is not sufficient to explain the 5(3) to 5(7) line that 
we observe. This emission requires higher gas temperatures that are reached only in faster 
shocks. Therefore we have tried to match the observed column densities with a combination 
of slow and fast shock components. Both the pre-shock density and initial o'pr were kept 
the same as for the slow shock component, and the velocity of the fast shock component was 
adjusted in order to match the 5(3) to 5(7) lines. We found that the emission of these lines 
were well matched for a shock velocity of 36-53 km/s and a filling factor of typically 1%. 
In fact, this two component model provides an excellent fit to the observations (see Fig. [s]), 
with the exception of the 5(0) line [E^p = 510 K] in IRAS-2, IRAS-7 and SVS13B-S, which 
is somewhat underestimated by the model. Possibly, this line arises in a separate, non-shock 



component towards these sources. Indeed, in their study of supernovae remnants, Neufeld 



et al. (2007) found that the 5(0) emission was more extended than other spectral lines they 
observed. In one source, they found evidence for both a shock excited 5(0) component and 
a diffuse emission component. 

Table |3] gives the best-fit model parameters for the slow and fast shock components. 
Towards all positions, we obtain a good fit for a pre-shock density of 10^ cm~^. The shock 
velocity in the slow shock component is between 12 and 24 km/s, while that in the fast shock 
component is between 36 and 53 km/s. For all sources but IRAS 4-SW, the emission is well 
matched assuming an initial o'pr of 0.01. For the latter, we obtain a slightly better fit for 
an initial o'pr of 1.0. However, the grid of models covers only a relatively small sample of 
o'pri (0.01, 1, 2 and 3) so this parameter is not constrained precisely. It is likely that oprj 
values ranging between 0.01 and 1 would also provide a good fit to the observations. The 
filling factor ranges from 2 to 15% for the slow component, and from 0.8 to 3% for the fast 
component. Since shocks driven by a coUimated jet are not planar but have a bow-shape, 
it seems natural that we observe a range of velocities: at the apex of the bow, the shock is 



moving faster with respect to the wings of the bow (Smith et al. 1991). This high- velocity 
component has a small beam filling factor, because the high velocities are present only at 
the apex of the bow. On the other hand, slower gas - leading to lower excitation - is present 
in the broader wings of the bow. 

It is worth noting that the fit is not be unique. For example, we found that the emission 
observed in HH 7 can also be fitted by a model with an higher pre-shock density (10^ 
cm~^), but a smaller pre-shock velocity (15 and 30 km/s for the slow and fast component. 
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respectively). In this case, the observations require an initial opr of 1, because the para to 
ortho conversion at these shock velocities is not efficient enough. Another reasonable fit was 
found for a pre-shock density of 10^ cm~^ and pre-shock velocities of 10 and 20 km/s for the 
slow and fast component, respectively. Therefore the pre-shock density and the pre-shock 
velocity are difficult to constrain simultaneously from the present observations. However, 
other constraints on these parameters have been obtained using FIR lines observations. 



Molinari et al. (2000) observed CO and H2 lines with ISO, and compared these observations 
with the predictions of Kaufman & Neufeld (1996) shock models. They found that the 



observed emission was well reproduced by a shock with a velocity ranging between 15 and 
20 km/s, and a pre-shock density of 10^ cm~^. Giannini et al. (2001) also compared the H2O 



and CO line emission observed with ISO-LWS towards the red-lobe of the IRAS 4 outflow to 



the model predictions of Kaufman & Neufeld (1996). Their observations are consistent with 
a velocity ranging between 15 and 20 km/s, and a pre-shock density of 10^ cm~^ (see their 
Fig. 8). Thus the pre-shock density adopted here is consistent with the values derived by 
these authors using CO emission, which is a much more sensitive probe of density. It is also 
interesting to note that the shock velocity range they obtain agree well with what we obtain 



here for the slow shock component. Furthermore Smith et al. (2003) obtained images of the 
2 — 1 S{1) and 1 — S{1) H2 emission lines around the HH 7 bow shock. They found that 
their observations were well reproduced by a C-type shock model with a pre-shock density 
of 8 X 10^ cm~^, a bow speed of 55 km/s and a magnetic field of 97 fiG (corresponding to 
6 ~ 1 for that pre-shock density). These parameters agree remarkably well with the values 
obtained here. The bow speed is slightly higher than the velocity we derive for the fast shock 
component (45 km/s). This is consistent with models that predict that pure rotational H2 
lines arise in the wings of the bow, where the shock velocity is smaller than at the apex (see 



Smith et al. 2003 Fig 5). 



An important result of the present study is that the H2 observations cannot be repro- 
duced by an initial ortho-to-para ratio greater than 1; models with opvi = 2 or 3 produce 
"fiat" rotational diagrams, i.e. opr values that are close to the high temperature limit of 
3. This indicates that H2 in the cold gas (i.e. before the passage of a shock) is mostly in 



para form. Flower et al. (2006) studied the evolution of the opr as a function of time in a 



pre-stellar core. In their computations, H2 is assumed to form on the grains with an opr of 
3. Reactive collisions between 0-H2 and p-Hj]" (or o-Hj]") forms P-H2, and consequently the 
opr decreases with time. The steady state value (~ 3 x 10"^ for T = 10 K, and ~ 3 x 10~^ 
for T = 30 K) is reached after t ~ 3 x 10^ yr, for a density of 10^ cm~^. Our observations 
and modelling suggest that the initial ortho-to-para ratio is lower than 1. According to the 



Flower et al. (2006) computations, an ortho-to-para of 1 is reached after 3 x 10 year (see 



their Fig. 1); this is comparable to the life-time of nearby molecular clouds obtained by 
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Hartmann et al. (2001). Therefore the value of the initial opr is consistent with that is 



expected from chemical models. 



5.2. Thermal history of the shocked gas 



The measure of the opr along the several outflows in the cloud provides constraints on 
the thermal history of the gas. Our observations and models suggest that in the cold gas, H2 
is mostly in para form. Shock waves heat up the gas and trigger para-to-ortho conversion by 
reactive collisions with H, so the opr increases. However, the reaction has a large activation 
barrier (~3900 K), so it is only efficient for gas temperatures greater than a few thousand 
Kelvins. As the post-shock gas cools down, the interconversion becomes less efficient, and 
the opr remains frozen. Of course, reactions with B.^ will tend to decrease the opr in the 
post-shock gas. However, this occurs on relatively long timescales (~ 10^ yr for nn = 10^; 



Flower et al. 2006), when compared to the dynamical timescale of the flows. Thus the opr 



retains a memory of the thermal history of the gas. 

Constraints on the thermal history of the gas may be placed by simple considerations. 
If we assume that a parcel of the gas is warmed-up by a passage of a shock wave to a 
temperature T for a duration r, and that ortho-para interconversion occurs only by reactive 



collisions with H, then the fraction of H2 in ortho form in this parcel is given by (Neufeld 



et al. 2006): 



/ortho (t) = /ortho exp 



-n(H) fcr 

^LTE 
J ortho 



/■LTE 
J ortho 



exp 



-n(H)fcr 

fLTE 
J ortho 



where /^tho the initial ortho fraction, /^™ is the ortho fraction at the local thermodynamic 

sum of the rate coefficients 



equilibrium, n(H) is the hydrogen density, and k 
for para-to-ortho and ortho-to-para conversion, given by k 



10~^^exp(-3900/T) cm^ s~^ 



^Schofleld 


1967 




Following 



10 



-3 



mainly produced by the formation of water from atomic oxygen, so that n(H)/n(H2 
Thus, for a pre-shock H2 density of lO'' cm~'^, n(H) ~ 10 cm~'^. Fig. 10 shows the opr as 
a function of time for different gas temperatures, assuming an initial opr of 0.01. In HH 7, 
we measure an opr of 0.37 and 1.99 in the warm (T = 611 K) and hot (T = 1401 K) 



^In the analogous equation given by (Neufeld et al 



2006), fcpo, the rate coefficient for para-to-ortho 



conversion alone, was erroneously given in place of k. The former is 25% smaller than the latter (a difference 
which is negligible relative to other uncertainties). 
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components, respectively. This corresponds to timescales of 
the warm and hot components, respectively. 



8 X 10^ and ~ xlO^ yr for 



The values we obtain can be compared to the prediction of shock models. The timescale 
during which the gas temperature is elevated on passing through a non dissociative shock 
is given by Ns/{nQVs), where Ng is the column density of the warm shocked gas, hq is the 



pre-shock density, and Vs is the pre-shock velocity (Neufeld et al. 2006). Assuming a pre- 
shock density of 10^ cm~^, and a pre-shock velocity of 20 km s~^, we obtain, using the 



analytical expression for Ns from Neufeld et al. (2006), a shock timescale of 400 yr. This is 



about an order of magniture lower than the timescales we obtain for both the warm and hot 
components of HH 7. This may suggest that several shocks waves have passed through the 
quiescent gas and have progressively increased the opr, so opvi is greater than 0.01. If we 
assume that opvi is 0.37 - the value we observe in the warm component of HH 7, at the tip 
of the HH 7-11 outflow - we found, using Eq. g, that 800 yr are needed to reach the value 
of 1.99 measured in the hot component, assuming a temperature of 1401 K. This is broadly 
consistent with the 400 yr timescale obtain above. 



5.3. Global outflow properties 

We have obtained 5.2 - 35.8 /im spectroscopic maps encompassing nearly the full spa- 
tial extent of 5 molecular outflows from embedded young protostars (HH 7-11, SVS 13B, 
IRAS 4A, IRAS 2, IRAS 7). In the following section we will explore how the information 
obtained in these maps provides a direct measurement of the mass loss rate and outflow 
energetics from these young stars. 



5.3.1. Determination of total outflow energy loss 

As shown in Fig. 4 a-i, we have detected optically thin emission from eight rotational 
transitions of H2 within these flows. Using these data, we can calculate the total H2 luminos- 
ity in rotational emission, Lhj, in the outflow by summing the emission over all rotational 
states and at each position within the region associated with each outflow. The transitions 
we have observed account for most of the luminosity carried by H2 rotational emissions, the 
brightest transitions falling within the range with Jup = 2 — 9 covered by Spitzer-IRS. We 
provide the value of for each flow in Table |4j The H2 luminosity is related to the total 
outflow cooling rate by Ltot = ^-^H2) where /c is the fraction of the total cooling due to H2 
rotational emissions. 
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The shock models of Kaufman & Neufeld ( 1996 ) can be used to explore the contribution 
of H2 to the overall cooling of the outflow via emission. Based on analysis of the H2 emission 
(§ 5.1 ) we find typical shock parameters of riQ ~ 10"^ cm~^ and Vg ~ 10 — 50 km/s. Under these 
conditions Kaufman & Neufeld (1996) find that fc ~ 0.2 — 0.7. The majority of H2 cooling 
is via rotational lines, with only few percent contribution from H2 vibrational emission. The 
rest of the cooling flux arises from rotational, and to a lesser extent vibrational, lines of H2O 
and CO, along with [O I] fine structure emission. 

The Infrared Space Observatory has observed the primary cooling lines of CO, [O I], 



H2O, OH, and H2 with the SWS and LWS spectrometers in some of our sources. Molinari 



et al. (2000) observed HH 7 and compile the luminosities of the major cooling lines within 
the bandpass. For the HH 7 shock they find that 4.8 x 10"^ Lq is released by O I, 2 x 10"^ Lq 
by CO, 0.7 X 10-2 Lq by H2O, and 2.3 x 10"^ Lq by H2. However, the O I, CO, and H2O 
observations were obtained using the LWS instrument with an 80" beam centered on HH 7, 
while the H2 observations covered only the S{1) — S{5) transitions within a 14" x 20" SWS 
beam. Our observations in Fig. 4a-i demonstrate that the H2 emission extends well beyond 
the SWS aperture and the ISO SWS value is therefore a lower limit. We have convolved our 
data with a simulated 80" LWS beam and have recomputed the total luminosity of H2 from 
5(0) — S{7). We find that the H2 luminosity in an 80" beam centered on HH 7 is 0.1 Lq. 
Note that in this source, the luminosity of the H2 1 — S{0) ro- vibration line is 3 x 10~^ Lq 
(Khanzadyan et al. 2003), well below the H2 luminosity in rotational emission. Combining 
the ISO observations with this revised value, we find that the H2 provides about 50% of the 
total cooling (/c ~ 0.5), in good agreement with models. 

Information for other flows can be gleaned from the summary of ISO cooling lines of 



Giannini et al. (2001) which included two of our sources: the IRAS 2 north-south flow, and 



the IRAS 4A outflow. For these sources, we find (after adapting the results of Giannini et al. 
to our adopted distance of 220 pc) that fc ~ 0.25. Using this information the total cooling 
luminosity of the outflows in NGC 1333 is approximately Ltot ~ 2 — 4 Lhj- For HH 7-11 we 
will adopt a fc = 0.5 and for all other flows (which are not as prominent as HH 7-11) we 
adopt fc = 0.25. 



5.3.2. Determination of the stellar mass loss by outflows 
The total energy loss generated by the outflow driven shock determined above is a 



direct measurement of the mechanical luminosity of the outflows. In this fashion, ^MqV 



2 



(1 — fm)Ltot = (1 ~ fm)Y-^^2y where (1 — fm) is the fraction of shock mechanical energy 
that is translated into internal excitation, as opposed to translational motion. [Kaufman &: 



Neufeld (1996) estimate that (1 — Z^) ~ 0.75 and we can then use the H2 coohng emission 
to determine the outflow mass loss rate, from the young star. This equation can be 
simplified and placed in terms of the observed H2 luminosity and estimated shock velocity: 



fc VlO-2 LqJ VlOkm/s 



M© yr-i. 



(2) 



To determine the shock velocity, we use the models of H2 emission discussed in § 5.1 In 



general, the data is best fit by a combination of fast and slow shocks with the slow shock 
encompassing a larger filling factor within the beam. We have therefore estimated the average 
shock velocity for each flow by using a filling factor weighted average of the shock velocity 
based upon the model fits. In Table |4] we provide the observed H2 luminosity, assumed shock 
velocity, and derived outflow properties. 

In general we find that the outflows from the Class protostars in the NGC 1333 cloud 
have values M^, ~ 0.6 — 2 x 10~^ Mq yr~^. These values are comparable to values estimated in 
the literature (e.g. Bontemps et al.|[l996 Giannini et al.|[2C)01 Hatchell et al.|[2007 ). In some 



cases measurements include sources within our sample. For example, Giannini et al. (2001) 
use FIR cooling lines towards IRAS 4 to derive = 0.4 — 1.4 x 10~^ Mq yr~^ where we 
estimate M^ = 2x 10"^ Mq yr'^ The derived mass outflow rates can be related to the mass 
accretion rate (Ma) using theoretical models of outflow generation. While the details and 
results can vary, a typical number is Mo ~ 0.1 Ma (Pudritz et al. 2007), and our estimates 
M0 yr^^. In this regard, di Francesco et al. (2001) independently derived 



provide Ma ~ 10 

a mass infall rate of 1 x 10"^ Mq yr"^ towards IRAS 4A using P Cygni spectral line profiles 
of H2CO. In addition, Maret et al. (2002) derived an accretion rate of 5 x 10~^ Mq yr~^ 
by modelling the FIR water lines emission from IRAS 4A and 4B envelopes. This is in 
reasonable agreement with our estimate of 2 x 10~^ Mq yr~^ for this source. 

It is useful to place our results in a broader context. Previous estimates of these quan- 
tities for molecular outflows have generally used velocity-resolved CO emission to provide 
a measure of the total H2 mass in the flow, and used the velocity extent of the lines to 



trace the kinematics (e.g. Lada 1985 Bontemps et al. 1996). However, there are a number 



of uncertainties inherent in this method. The CO abundance has been calibrated to H2 in 



some clouds (e.g. Frerking et al. 1982), but can vary if some CO remains frozen on grains 



e.g. Bergin & Tafalla 2007) or is photodissociated by radiation generated in the shock or 



leaking through the outflow cavity. In addition, the CO molecules can be excited in both 
the low temperature natal core and in the warmer outflow. Thus towards the central core 
the spectral line contains a mixture of emission from the high velocity dispersion core and 
low dispersion core. To some extent these components can be disentangled, but there is 
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significant uncertainty regarding the outflow mass present at low velocities that also corre- 
spond to the core dLadafc Fich]p96| . CO emission itself is often optically thick, even in the 
extended line wings. Observations of isotopologues are therefore required, but not always 



available to correct for the effects of optical depth (e.g. Yu et al. 1999 Arce &; Goodman 



2001). Finally, estimates of the CO velocities are dependent on the outflow inclination. In a 



handful of cases, this can be directly measured if proper motion data can be used to measure 
the tangential velocity to compare to the observed radial velocity in the spectral lines. In 
the majority of cases this information is not available, and inclinations have instead been 
estimated by models of the flow velocity at different inclinations compared to the observed 



emission distribution (e.g. Cabrit et al. 1988). However, this estimate is also uncertain. 



Our measurements do not suffer from any of these uncertainties. The H2 lines are 
optically thin and our observations encompass the most emissive II2 rotational lines, which 
we simply sum over rotational states and positions to derive the total emission. We have 
used shock models to calibrate the amount of cooling accounted for by other molecules. In 3 
sources covered by our maps, we have conflrmed that the models are accurately reproducing 
the observed distribution of cooling power. Since H2 is the dominant species and a primary 
coolant we are directly tracing the energy loss in the most abundant species. H2 does not 
emit within the core itself and thus we have no line of sight contamination. By using the 
cooling luminosity to trace outflow energetics our measurements are completely independent 
of the outflow inclination. As an example our momentum flux estimates are a factor of 
~ 4 — 8 above values - uncorrected for inclination - derived from CO emission in same flows 
in NGC 1333 ( Hatchell et al.||2007 ). On the other hand, our estimates are in good agreement 



the values (corrected for inclination) of Sandell & Knee (2001). 



Our conclusion is similar to that of Giannini et al. (2001) who used ISO observations 



of FIR cooling lines as a measure of the total mechanical luminosity. The major difference 
is that we use spatially resolved observations of H2 and supplement our results with those 



of Giannini et al. (2001) to derive the total cooling loss. Given that we have observed the 



major coolants, the primary uncertainty in the mass loss rate derivation is the assumed shock 
velocity. The expected range of shock strengths to create the observed H2 emission is ~10 - 
50 km/s, with adopted values of ~ 20 — 30 km/s (Table |3]). We therefore estimate that the 
mass outflow rate values are accurate to within a factor of 3 (given a distance of 220 pc). 



5.3.3. Impact of flow on natal core 



It is well known that outflows inject signiflcant energy and momentum in the the sur- 
rounding cloud (e.g. Lada 1985 Arce et al. 2007; Davis et al. 2008) with some suggestions 
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that these flows drive supersonic motions ( Mac Low &: Klessen||2004 ) . It is also possible that 



the outflow is the key player in the disruption of the natal core (e.g. Tafalla & Myers 1997 



Fuente et al. 2002). Our data offer a new opportunity to explore the question of the outflow 



impact on the surrounding material in the specific case of low mass star formation. 

The total momentum injected by outflows into a core is given by P = P r^yn ~ 0.3 Mq 
km s~^, where P = Vg and r^j,„ is the outflow dynamical timescale. Using T^yn values from 
Knee & Sandell (2000), we find P values ranging from 0.1 to 0.4 Mq km s~^ (see Table [4|. 
If a similar level of outflow activity persists during the lifetime of the embedded phase of 
~ 5 X 10^ yrs (Evans et al. 2008) then the total momentum over this phase is ~ 4 — 20 
Mq km s~^. There is some evidence in the literature of a decline in the outflow momentum 



flux during the transition from Class to Class I (Bontemps et al. 1996) However, Hatchell 



& Fuller (2008), with a homogeneous sample of sources in Perseus, do not confirm this result 



finding similar momentum fluxes for sources with comparable luminosities and masses. 

Based on these estimates when the entrained outflowing material slows down to 1 km 
s~^ then it will have swept up 4 — 20 Mq of material. Assuming a typical core mass of ~ 5 Mq 
then the outflow is clearly capable of destroying the core during the lifetime of the embedded 
phase. Of course, the flow will only disperse material in its path and the typical outflow cone 
angle in these sources is ~ 60° (J0rgensen et al. 2007). However, observations suggest that 



the outflow cone angle increases with age (Velusamy & Langer 1998; Arce & Sargent 2006). 



Estimates of core masses vary by nearly an order of magnitude in the literature (Lefloch 



Arce & Sargent 


2006 


)• 


the literature ( 


L/efloch 



et al. 1998 Sandell & Knee 2001 Walsh et al. 2007) so deflnitive comparisons within our 



small sample are difficult. Nonetheless these data clearly suggest that the outflow could be 
the main mechanism for core dispersal. 

To support our assertion that the transfer of outflow momentum to the core is a main 



mechanism for core dispersal we can also compare the energetics. Tafalla & Myers (1997) 
demonstrated that the outflow from a massive star (Mon R2) exceeds or is comparable to 
the gravitational binding energy of its core. Thus the Mon R2 outflow has likely sculpted 
that object. However massive star outflows are more energetic than the low mass objects 
sampled in our maps and this result may not scale to low mass systems. The average 
total energy generated by the current generation of outflows in the center of NGC 1333 is 
{Lii,J fc X Tdyr) = 2 X 10^^ ergs. The gravitational binding energy of a sphere is aGM'^/R 
with a = 1 for cores with a density proflle (MacLaren et al. 1988), as expected for 
these embedded protostars. For a typical core mass of ~ 5 Mq and radius of 0.03 pc the 
binding energy is ~ 4 x 10^^ ergs. If we assume that a similar level of activity exists over the 
lifetime of the embedded phase then the outflow energetics are equal to or even exceed the 
gravitational energy of the core. In sum, both the flow energetics and momentum suggest 



- 20 - 



that outflows are the primary method for core dispersal. Molecular gas observations show 
that NGC 1333 cores discussed here are embedded in a much larger but lower density cloud 
(e.g. Pineda et al. 2008). As typical for GMCs, the cloud binding energy is much greater 
than those of each individual cores. Thus the flows which disrupt the cores are not major 
players in cloud disruption as the cloud momentum and energy far exceed that provided by 
the current generation of flows over their lifetime. 



6. Conclusions 

We have presented Spitzer-IRS maps of seven pure rotational H2 lines in the NGC 1333 
star forming region. These observations cover a region roughly 20' x 13' that encompass a 
number of YSOs and outflows. We have analysed these observations using the rotational 
diagram technique in order to derive the excitation temperature and the H2 ortho-to-para 
ratio. Furthermore, these observations have been compared with the predictions of a planar 
stationary shock model. Finally, we have used the H2 total luminosity to estimate the mass 
outflow rates, the mass infall rate onto the central objects, and the momentum and energy 
injected into the cores and the cloud. Our main conclusions are as follow. 

• H2 line emission is detected along several outflows in NGC 1333. In particular, we 
detect line emission along both lobes of the HH 7-11 outflow, along the south-west 
lobe of the outflow driven by IRAS 4A, along the north, south and east lobes of the 
IRAS 2 outflow, as well as around IRAS 7. In addition, we detect some faint emission 
south of SVS 13, which is associated with the outflow originating from SVS 13B. Diffuse 
emission is detected north of this source, and might be associated with this outflow as 
well. 

• A rotational diagram analysis on several aperture averaged spectra indicate the pres- 
ence of at least two gas components towards each position. The warm component has 
a rotational temperature between ~ 300 and 600 K and an opr between ~ 0.3 and 
0.7, while the hot component has a rotational temperature between ~ 1000 and 1500 
K and an opr between ~ 1.9 and 2.2. 

• Comparison of the line fluxes measured towards these positions with the predictions of a 
planar stationary shock model indicate the presence of slow and fast shock components. 
A good flt of the observation is obtained for an initial opr of 0.01 (except in one source 
where the data are better flt with an initial opr of 1), a pre-shock density of 10^, 
and a shock velocity of ~20 km/s and ~45 km/s the slow and fast shock components, 
respectively. This indicates, in agreement with earlier studies, that H2 is mostly in para 
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form in dense molecular clouds (i.e. opr < 1). This is consistent with the predictions 
of chemical models, provided that the cloud lifetime exeed a few million years. 

• Maps of the rotational temperature and the opr in warm components show some 
important variations of these quantities across the mapped region. We confirm the 
presence of gas with a relatively high rotational temperature and low opr at the tip 
of the HH 7-11 outflow. This gas is probably "fresh gas" with a low initial opr that 
has been heated-up by the passage of a shock, but whose opr hasn't had time to 
increase significantly. The same feature is observed on both lobes of the IRAS 7 flow, 
although less clearly. Spatial averages of the opr and the rotational temperature along 
the flow support this hypothesis. Interestingly, no spatial correlation between the two 
quantities is found. 

• Comparison between the total H2 luminosity and the luminosity of other lines observed 
by ISO in HH 7-11 and IRAS 4 indicates that between a fourth and half of the total 
cooling in these outflows occurs through the H2 lines. Therefore H2 lines can be 
used to estimate the kinetic energy injected into the natal cloud by these flows, as 
well as outflow mass loss rate which are found to be 0.6 — 2 x 10~^ Mq yr~^. The 
latter quantity also places indirect constraints on the mass accretion rate onto the 
protostar: dynamical models predict that the outflow mass loss rate is typically 10 
times lower than the accretion rate. Using the dynamical timescale obtained from 
millimeter observations, we derived the total momentum injected in the clouds by the 
outflows, and compared it to the core and cloud binding energy. We found that outflows 
have the potential to disrupt individual cores, but are probably not major players in 
cloud disruption. 



This work is based on observations made with the Spitzer Space Telescope, which is 
operated by the Jet Propulsion Laboratory, California Institute of Technology under a con- 
tract with NASA. Support for this work was provided by NASA through an award issued 
by JPL/Caltech for program 7^20378. 

Facilities: Spitzer (IRS) 

A. Variation of the opr with the temperature 



As mentioned in § 4^, some important spatial variations in both the opr and the 
rotational temperature are observed in NGC 1333. In particular, we observe, at a tip of 
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the HH 7-11 outflow, a bow shock-like region with a low o'pr (~ 0.5) and high rotational 



temperature (~ 650 K; see Fig. 11). This low o'pr region wraps nicely around the H2 5'(1) 
rotational line emission, and suggests that it corresponds to gas which has been recently 
heated-up by the passage of a shock, but whose o'pr has not had time to reach the equilibrium 
value. Interestingly, this region also corresponds to the part of the flow where the rotational 
temperature is the highest. Around IRAS 7, we also observe low o'pr regions around the H2 



5'(1) emission contours (see Fig. 12), but variations from one pixel to the other are important 



and it is less clear if this region correspond to high temperature gas, like in HH 7-11. 

It is interesting to compare the variations of the opr and the rotational temperature as 
a function of the position of the flow in a more quantitative fashion. In order to do this, we 
have computed the average of the opr and rotational temperature along the outflows that 
show some degree of symmetry, and where both quantities were measured in a sufficient 
number of pixels. In order to "capture" the spatial variations, we have defined a series of 
ellipses of the same eccentricity, oriented along the fiow and with an apex that is centered 
on the outfiow driving source. The eccentricity of the ellipses has been chosen in order to 
match the shape of each outfiow lobe as observed in the H2 5(1) emission. We have then 
computed the average of the opr ratio and the rotational temperature in the region between 



two consecutive ellipses (that is two ellipses with different major axes). On Fig. 11 and[l2| 
we show how those ellipses where defined for HH 7-11, IRAS 7 (east and west lobes) and 
IRAS 2 (south lobe). We have not computed the average for the other outfiow lobes, because 
the opr ratio and the rotational temperature had greater uncertainties. 



Fig. 13 shows the average opr and rotational temperature, as a function of the ellipse 
semi-major axis. Along the HH 7-11 outfiow, we see a sharp drop in the average opr ratio 
35" from SVS 13. Interestingly, the temperature does not seem to be correlated with the 
opr. It increases slightly and 35" (although at a Icr level), and then starts to decrease. 
The highest temperature is reached where the opr ratio is minimal. Presumably, the hottest 
point corresponds to the shock front. In the post-shock gas, the para to ortho conversion 
has increased the opr^ but the conversion is incomplete and the ratio does not reach the 
LTE value. In the pre-shock gas the opr has also probably increased with respect to its 
initial value, but has not had time to reach the same value that we see in the post-shock 
gas (i.e. ~2). Note that for a pre-shock density of 10^ cm~^, and a shock velocity of 20 
km/s, the model predicts a shock front width - defined as the region of the shock where the 
temperature is greater than 1000 - of ~ 10'^ AU, i.e. about 5" at the distance of NGC 1333. 

The interpretation of the opr and rotational temperature along the other fiows is less 
straightforward. We see an increase of the rotational temperature along the IRAS 4 lobe, but 
the maximum temperature is reached at the tip of the outfiow; farther away no H2 emission 
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is detected and the rotational temperature and the opr could not be determined. Unlike 
the HH 7-11, the average opr ratio is roughly constant, and is ~ 2. The temperature in the 
east lobe of IRAS 7 flow seems to slightly decrease (note that the flrst ellipse contains only 
a couple of points, so the average is not very meaningful in this bin), while the opr is almost 
constant. On the other hand, the opr ratio and temperature may vary in a similar fashion 
that in HH 7-11. The opr drops sharply between 5 and 15", while the temperature seems to 
slightly increase in the same region. However, further away the temperature appears to be 
roughly constant. 
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6 8 10 12 14 



X (/xm) 

Fig. 1. — Averaged spectra observed with the SL module for 15" Gaussian apertures 
centered, from top to bottom, on HH 7 (a = 03^29"^08.45% 6 = +3ri5'29.2"; J2000) 
IRAS 2A {a = 03^28'^55.59^ 6 = +31°14'37.3"; J2000), IRAS 2A-S {a = 03^28"^54.04^ 
6 = +31°13'31.2"; J2000), IRAS 4A {a = 03''29°'10.29^ 5 = +31°13'31.8"; J2000), 
IRAS 4A-SW {a = 03''29°^06.62^ 6 = +31°12'17.7"; J2000), IRAS 4B {a = 03^29'^11.99^ 
6 = +31°13'08.9"; J2000), IRAS 7 (a = 03^29'^11.3P, S = +31°18'31.1"; J2000), 
IRAS 8 (a = 03^29'"12.5^ 5 = +31°22'08"; J2000), and SVS 13B-S {a = 03'^29'"05.0^ 
S = +31°13'14"; J2000). Pure rotational H2 transitions are indicated, as well as the 6.2, 7.2, 
8.6, 11.2 and 12.8 /im PAH features, and the 9.7 fim silicates absorption band. The spectra 
between ~ 5 and 7.5 fim and between ~ 7.5 and 14 /xm correspond to the second and first 
order of the module, respectively. 
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Fig. 2. — Same as Fig. [T]for the SH module. The position of the 15.2 /i CO2 ice feature is 
indicated. For clarity the spectral resolution has been degraded by a factor 4. 
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Fig. 3. — Same as Fig. [T]for the LH module. For clarity the spectral resolution has been 
degraded by a factor 4. The steps seen in some of these spectra are caused by offsets in 
response of the different orders of the module. 
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Fig. 4 a. — H2 5'(1) continuum subtracted map obtained with the SH module. YSOs are 
indicated by open stars, while several HH objects are indicated by open squares. The outflows 
discussed in the text are indicated by dashed arrows. The two positions along the IRAS 2 
and IRAS 4A outflows that are discussed in the text, IRAS 2A-N and IRAS 4A-SW, are 
also shown. Contours show the 3, 5 and 10 cr noise levels. See the electronic edition of the 
Journal for panels d-i. 
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Fig. 4 b. — H2 5(2) map obtained with the SL module. Contours show the 3 and 5 a noise 
level. 
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Fig. 4 c. — H2 5(3) map obtained with the SL module. Contours show the 3 and 5 a noise 
level. 
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Fig. 5. — H2 rotational diagrams obtained for 15" Gaussian apertures centered on HH 7, 
IRAS 4A-SW, IRAS 2A-S, IRAS 7 and SVS 13B-S. The black open circles correspond to 
the observations. The black open triangle is a la upper limit. The black solid lines are 
rotational diagram obtained assuming that two gas components, a warm and a hot one, are 



present (see § 4.1). 
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Table 2. Rotational diagrams results 



Warm component Hot component 



Source 


iV(H2) 


^rot 


opr 


iV(H2) 


^rot 


opr 




(cm~^) 


(K) 




(cm~^) 


(K) 




HH 7 


5.9 X 10^9 


611 


0.37 


6.2 X 10^^ 


1401 


1.99 


IRAS 4A-SW 


5.9 X 10^9 


342 


0.69 


3.4 X lO^s 


1513 


1.94 


IRAS 2A-N 


2.0 X 10^9 


371 


0.33 


5.6 X 10^^ 


1055 


2.09 


IRAS 7 


4.2 X 10^^ 


300 


0.31 


3.3 X 10^^ 


1268 


2.11 


SVS 13B-S 


3.2 X 10^9 


312 


0.51 


1.5 X lO^s 


1337 


1.43 



Table 3. Best fit shock model parameters for aperture-averaged positions 





Slow shock component 




Fast shock component 


Source 


n(H) 


V 


opn 


// 


n(H) 


V 


opvi 


.// 




(cm~^) 


(km s~^) 




(%) 


(cm~^) 


(km s~^) 




(%) 
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10^ 


22 


0.01 


15 


10^ 


45 


0.01 


3 
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10^ 
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1 


2 


IRAS 2A-S 


10^ 


23 


0.01 


4 


10^ 


50 


0.01 


1.5 


IRAS 7 
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0.01 


2.2 


10^ 


53 


0.01 


1.1 


SVS 13B-S 


10^ 


14 


0.01 


3 


10^ 


36 


0.01 


0.8 
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Fig. 6. — H2 column density map (color points) on top of the H2 S{1) emission (black 
contours). Regions of the map that appear in white correspond to points where a rotational 
diagram could not be constructed because of missing data or insufficient signal-to-noise ratio. 



Contour levels are those of the 5'(1) emission taken from Fig. 4 a . Pixel size is 5". 
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Fig. 7. — Same as Fig. |6]for the H2 rotational temperature. 
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Fig. 8. — Same as Fig. |6]for the H2 ortho-to-para ratio. 
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Fig. 9. — Same as in Fig. |5} but with the two velocity components shock model predictions 



over plotted (see § 5.1 ). The blue and red squares connected by dotted lines show the emission 
from the slow and fast shock components, respectively. The orange squares connected by 
dotted lines show the total emission from these two components. Note that some of the blue 
squares (at low energies), as well as some of the red squares (at high energies) are masked 
by the orange squares. 
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Fig. 10. — H2 ortho-to-para ratio as a function of time for different gas temperatures. An 
initial opr of 0.01 and H density of 10 cm~'^ are assumed for all curves. 



Table 4. Outflow properties 



Source 


Lh2 


fc 


Vs 






P 


Tdyn 


P 




(10-2 Lo) 




(km s-"*^) 


(Mo yr- 




(Mq yr-1 km s-^) 


(yr) 


{Mq km s-^) 


HH 7-11 


11.1 


0.50 


26 


2 xlO^ 


-6 


5 xlO^^ 
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0.29 


IRAS 4A 
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0.25 


16 


2 xlO" 


'6 
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0.34 


IRAS 7 


1.8 


0.25 


34 


5 xlO" 


-7 
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0.11 
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0.8 


0.25 
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^From [Knee fc San"deII| ( [2000| 

^North-south outflow only 
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Fig. 11. — Side-by-side comparison between the rotational temperature and opr around 
SVS 13 (upper panels) and IRAS 4 (lower panels). The position of several HH objects 



around SVS 13, from Walawender et al. (2005), are indicated . The ellipses that were used 



to average the opr and the rotational temperature (see § \M are shown. 
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Fig. 12. — Same as Fig. 



11 



for the IRAS 2 (upper panels) and IRAS 6 (lower panels). 
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Fig. 13. — Average H2 ortho-to-para ratio (blue histograms in solid lines) and rotational 
temperature (red histograms) between consecutive ellipses of the same eccentricity, as a 
function of the ellipse semi-major axis. The error bars are 1 a standard deviations computed 
from the variance in each bin (i.e. the region between two ellipses). The blue histogram in 
dotted lines show the maximum and minimum values of the opr in each bin. 
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Fig. 4 d. — H2 5'(0) map obtained with the LH module. Contour levels are drawn at 3 cr. 
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Fig. 4 e. — H2 S'(2) map obtained with the SH module. Contours show the 3 and 5 a noise 
levels. 
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Fig. 4 f. — H2 5'(4) map obtained with the SL module. Contours show the 3 a noise level. 



- 48 - 



31°20'00" 



ario'oo" 



NGC 1333 Hs S(5) 




f 



svs 



HH 



mm- - svs 



IRAS 2A' 



IHAS 2B 



-SVS 13B--S1 



. ^^i^IRAS''4^-SW 




3^29^00^ 



10 



-3 



CD 

P 
M 



I 

to 



-4 ^ 

10 I 



a (J2000) 

Fig. 4 g. — H2 S{5) map obtained with the SL module. Contours show the 3 cr noise level. 
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Fig. 4 h. — H2 S{6) map obtained with the SL module. Contours show the 3 a noise level. 
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Fig. 4 i. — 



H2 S{7) map obtained with the SL module. Contours show the 3 a noise leveL 



